Initialisation

## [1] "number of cores available = 1"
#Phi[1] ; eta = valeur de fin  Phi[2] = valeur du noeud  Phi[3] = echelle
m <- function(t, eta, phi) (phi[,1] + eta)/(1+exp((phi[,2]-t)/phi[,3]))
#=======================================#
param <- list(sigma2 = 0.05,
              rho2 = 0.1,
              mu = c(5,4,0.5),
              omega2 = c(0.5,0.1,0.01))
#=======================================#
t <- seq(2,6, length.out = 10) #value of times

# --- longitudinal data --- #
dt_NLME <- NLME_data(G = 40, ng = 12, time = t, fct = m, param = param)

getDim(dt_NLME)
##    G   ng    n    N   F. 
##   40   12 4800  480    3
Y <- dt_NLME$obs

SAEM avec simulation par MCMC

\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]

Phi <- fct_vector(function(sigma2, rho2, mu, omega2) 
  c(-n/(2*sigma2), -N/(2*rho2), -G/(2*omega2), G*mu/omega2), dim = c(1,1,F.,F.) )$eval

S <- fct_vector(function(eta, phi) mean((Y - get_obs(m, dt_NLME, eta = eta, phi = phi) )^2 ),
              function(eta, phi) mean(eta^2),
              function(eta, phi) apply(phi^2, 2, mean),
              function(eta, phi) apply(phi  , 2, mean), dim = c(1,1,F.,F.) )

Metropolis Hastings

loglik.phi <- function(x, eta, Phi)
{
  id <- c(1,3,4) #indice de S_1 et S_{3,.} puis S_{4,.}
  sum( Phi%a%id * S$eval(eta, x, i = id) )
}
loglik.eta <- function(x, phi, Phi) 
{ 
  id <- c(1,2)
  sum( Phi%a%id * S$eval(x, phi, i = id) )
}

SAEM

Initialisation

# ---  Initialisation des paramètres --- #
para <- param %>% lapply(function(x) x* runif(1, 1,1.2))
para$rho2 = 0.1 ; para$mu <- c(4,5,2) ; para$omega2 <- rep(.1,3)

# --- Initialisation des chaines MC : Z_0 ---
Z <- list(eta = list( rnorm(G*ng, 0, para$rho2)  %>% matrix(ncol = 1) ),
          phi = list( matrix(rnorm(F.*G, para$mu, para$omega2), nrow = F.) %>% t ) )

Etape simulation et maximisation du SEAM

sim <- function(niter, h, Phih, eta, phi)
{
  M <- length(phi)

  eta <- 1:M %>% lapply( function(i)
    MH_High_Dim_para_future(niter, eta[[i]], sd = .036, loglik.eta, phi[[i]], Phih, cores = 1))
  
  phi <- 1:M %>% lapply( function(i)
    MH_High_Dim_para_future(niter, phi[[i]], sd = c(.028, .036, .021), loglik.phi, eta[[i]], Phih, cores = 1 ))
  
  list(eta = eta, phi = phi)
}

maxi <- function(S)
{
  list(sigma2 = S%a%1,
       rho2 =   S%a%2,
       mu =     S%a%4,
       omega2 = S%a%3 - (S%a%4)^2 )
}
sigma2 rho2 mu1 mu2 mu3 omega21 omega22 omega23
Oracle 0.0514034 0.0973056 4.897666 4.002448 0.4832633 0.4364442 0.0730188 0.009188
Initialisation 0.0588703 0.1000000 4.000000 5.000000 2.0000000 0.1000000 0.1000000 0.100000 
niter <- 100
MH.iter <- function(k) ifelse(k<50, 20, 10)

res <- SAEM(niter, MH.iter, para, Phi, S$eval, Z, sim, maxi, burnin = 150, coef.burnin = 3/4, eps = 1e-3, verbatim = 2)
saveRDS(res, 'saem.rds')
gg <- plot(res, true.value = param)
## [1] "SAEM execution time = 11min 23sec"
Result of the SAEM-MCMC
sigma2 rho2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3
Real value 0.0500 0.1000 5.0000 4.0000 0.5000 0.5000 0.1000 0.0100
Estimated value 0.0506 0.1695 4.6820 3.9801 0.4839 0.3850 0.0747 0.0097
Rrmse 0.0125 0.6951 0.0636 0.0050 0.0323 0.2301 0.2527 0.0288